Clustering of inertial particles in turbulent flows 
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We consider inertial particles suspended in an incompressible turbulent flow. Due to inertia 
of particles, their velocity field acquires small compressible component. Its presence leads to a 
new qualitative effect — possibility of clustering. We show that this effect is significant for heavy 
particles, leading to strong fluctuations of the concentration. 
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I. INTRODUCTION 

Observing air bubbles in water or dust in air, one read- 
ily notices that inertial particles suspended in an inhomo- 
geneous flow tend to cluster. For example, such cluster- 
ing is widely used for flow visualization. Here we develop 
a statistical theory of this phenomenon. We describe the 
initial growth of concentration fluctuations and the sat- 
uration of that growth due to finite-size effects, imposed 
either by the Brownian motion or finite distance between 
the particles. Such theory is supposed to have numerous 
geophysical and astrophysical applications. A proper ac- 
count of concentration fluctuations is also necessary for 
a consistent theory of turbulent suspensions. 

Macroscopic description of a dilute suspension can be 
deduced from the the behavior of a single particle. Con- 
sider a small spherical particle with the radius a and the 
material density po suspended in a fluid with the density 
p. The particle's velocity v is related to the fluid velocity 
u by the equation dv/dt — (3du/dt — (u — v)/t s , where 
(3 = 3p/(p + 2p ) and r s = a 2 /(3^/3) is the Stokes time 
0. Both v and u are evaluated on the particle's trajec- 
tory q(t,r) that satisfies dtq = v and q(0,r) = r. The 
flow surrounding the particle is assumed to be viscous, 
which requires a r v , where r v is the viscous scale of 
the flow. This allows one to solve the system for v and 
q perturbatively in t s 



v = u + ((3- 1)t s [d t u + (u ■ V)u] . 



(1.1) 



If such particles are spatially distributed, it is possible to 
define the particles' velocity field v(t,r), which is com- 
pressible even if the fluid flow is incompressible [jl],||: 
(V • v) = (/3 - 1)t s V[(u • V)u]. Thus in the above ex- 
pansion we keep the terms up to the first term with non- 
vanishing divergence. As we show the smallness of this 
term may be compensated by large parameters, so that 
it can lead to significant effects. 

The concentration of the particles satisfies the 
diffusion-advection equation 



dtn + V(tm) = KV 2 n . 



(1.2) 



Every particle produces a relative perturbation of the 
flow that decays as an inverse distance from the parti- 
cle, i.e. as a/r. Since particles move coherently within 



the viscous scale r v , the condition a f v n(r)r~ 1 d 3 r ~ 
nar 2 , <C 1 has to be satisfied in order to be able to neglect 
their interaction. This condition is more restrictive than 



that of a small concentration, 



< 1. If 



< 1, the 



concentration field can be considered as passive, i.e. v is 
independent of n in Eq. (1.2). 

In statistically steady flows, velocity v in Eq. (1.2) 



must be considered as a random field with a station- 
ary statistics. Evolution of an arbitrary initial condition 
n(0,r) according to Eq. (1.2) ultimately results in the 
steady state of the concentration fluctuations. Making 
the decomposition n(0,r) = no + <5n(r), where no is the 
spatial average of n(0,r), one can write the solution in 
the form 



Jdr' G(t,r,r') + j ' dr' G(t,r,r')Sn(r') (1.3) 



n(t,r)=n 



We note that the Green's function, G, is nonnegative. 
At k = it is concentrated on the Lagrangian trajectory 
that passes through the observation point, r, at time t. 
At non-zero k, the Green's function is non-zero in some 
region around the Lagrangian trajectory. The size of that 
support region grows in time due to the combined ac- 
tion of the velocity and diffusion. As long as that size is 
smaller than the correlation length of the concentration, 
5n can be taken out of the integral. Then the expectation 
value of n depends on the ratio of no and the strength 
of initial fluctuations Sn. At large times when the sup- 
port of the Green's function becomes much larger than 
the correlation length, the second term contains contri- 
butions that cancel each other. The expectation value 
of n is then determined by the first term. Even though 
the second term in Eq. (1.3) may grow at these times, it 
is much smaller than the first term and only gives sub- 
leading dependencies. Thus the initial inhomogeneities of 
the concentration field arc irrelevant in studying the long- 
time evolution and the steady state. Therefore, we will 
consider a uniform initial concentration below. We will 
show that evolution of the uniform initial concentration 
distribution exemplifies most strikingly the inadequacy 
of the mean field picture for the problem. We choose the 
units so that no = 1. 

The paper is organized as follows. First, we analyze 
the initial period of growth when one can set k — (to 
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be referred to as the ideal case). We show that at this 
stage of evolution the moments of the concentration grow 
exponentially for quite an arbitrary velocity statistics (in 
case of time-decorrelated velocity this has been shown in 
H ) . The next section is devoted to the analysis of larger 
times. We show that diffusion modifies the growth of the 
moments and finally brings about saturation. 

II. THE IDEAL CASE 

The diffusivity of macroscopic particles is usually very 
small so that we will be interested in the statistics of n 
in the limit of small k. Starting with a uniform distri- 
bution one expects that at moderate times the diffusion 
term can be neglected until very thin structures are de- 
veloped. Let us consider this period of evolution. 

To find the concentration n{t\,rx) one has to count 
all the particles that come to a small volume (still, con- 
taining many particles) around r = r% and divide the 
result by this volume. To find which particles comes 
to this volume one has to track the trajectories of the 
particles backwards in time to t = 0. Particles per- 
form combined Lagrangian and Brownian motions. At 
t w ti, the size I of the region occupied by particles 
grows as l 2 {t) ~ K{t% — t). Velocity gradient A produces 
another mechanism of size stretching which takes over 
when / > \J k/X. Smallness of diffusion means that the 
Schmidt number, Sc = v / k is large so that \J k/X <C r v . 
As one goes further backwards in time the impact of dif- 
fusion on particles' motion becomes negligible. The time 
needed for I to reach the diffusive scale k/X is of the 
order of 1/A. On a larger time-scale the role of diffu- 
sion is to create an initial volume of finite size \J k/X. 
In other words, diffusion introduces the smallest scale 
into the problem, so that fluid particles cannot be lo- 
calized on the distances smaller than this scale. To find 
the concentration at a point in the ideal case one should 
track back an infinitesimal volume around that point. 
However, in the case of a nonzero diffusivity one should 
take the volume with a finite size \/ k/X. Since \J k/X 
can be smaller than the minimal scale of coarse-graining 
n^ 1 / 3 we introduce r^, which is given by the largest of 
two scales: J k/X and n -1 / 3 . To summarize, n(ii,ri) 
is given by the relative change of the volume of the size 
rd taken around the point r\ and tracked backwards in 
time along the Lagrangian trajectory. Let us note that 
the velocity gradient fluctuates in a random flow and so 
does the diffusion scale. We neglect these fluctuations 
(c.f. [Q) because they do not change the dependence of 
the concentration on large parameters, which arc either 
time in the transient regime or Reynolds and Schmidt 
numbers in the steady state. 

At t < A -1 ln(r v /rd) the rd-volume (the region that 
acquires the size ra at t = ti) always stays within the 
viscous scale and hence evolves in a uniform velocity 



gradient. Its relative change is the same as for an in- 
finitesimal volume, which means that the concentration 
behaves the same as in the ideal case. Equation (1.2) 



written in the Lagrangian frame therefore becomes the 
ordinary differential equation dn/dt = — n(V • v). Here 
(V • v) is a function of time, which fluctuates in a ran- 
dom flow. If the Lagrangian correlation time of the fluid 
velocity, it, is finite, which is true for most flows of geo- 
physical interest, then (V • v) ex V[(ix • V)u] has also a 
finite correlation time, r. At t 3> r, the concentration 
logarithm, X(t) = In [n(t)/n{Q)\ = -/ *(V • v)dt' , is a 
sum of a large number of random variables. The theory 
of large deviations assures that the probability density 
function (PDF) has the form V(X) oc exp\-ts(X/t)}, 
where s is a non- negative convex function ||. To cal- 
culate the moments of the concentration in the Eulerian 
frame one has to take every Lagrangian element with its 
own weight proportional t o its volume , i.e. to the inverse 
concentration (see Eqs. ( 2.20 ) and ( |2.2l| ) below). We 
obtain 

(n a (t,r))oc J dX exp[{a - 1)X -ts{X/t)\. (2.1) 

At large times, this integral can be found using the 
saddle-point approximation. The saddle-point X a is 
given by s'(X a /t) = a — 1, which implies X a cx t. 
Hence the moments generally behave exponentially in 
time: (n a (t)) oc exp(~j(a)t). 

Let us show that the conclusion on exponential behav- 
ior of moments is enough to establish the most interest- 
ing properties of this stage of evolution. The number 
of particles is conserved, i.e. (n) is time- independent. 
Hence 7(1) = 0. It is also obvious that 7(0) = 1. Due 
to the Holder inequality, the function 7(a) is convex. 
Therefore, 7(a) is positive for < a < 1 and nega- 
tive otherwise. Low-order moments decay whereas high- 
order and negative moments grow. The decay rate is 
(\og\n\)/t = ~dj(a)/da\ a= o < 0, i.e. n decays almost 
everywhere. Since the mean concentration is conserved, 
n has to grow in some (smaller and smaller) regions, 
which implies growth of high moments. The growth 
of passive scalar fluctuations in the particular case of a 
short-correlated compressible flow was described in || . 

The finite value of the coarse-grained volume comes 
into play at t > A -1 ln(r v /rd). Indeed, since particles 
separates backwards in time, the size of the volume ex- 
ceeds r v at t = 0. In other words, spots originated from 
different viscous domains come into contact at t = t\. 
A careful analysis of the time-span of the ideal case ap- 
proximation demands more detailed information on the 
Lagrangian dynamics at scales smaller than r v (later re- 
ferred to as small-scale Lagrangian dynamics). It is dis- 
cussed in subsection [I A. 
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A. Small-scale Lagrangian dynamics of compressible 
fluids with finite Lagrangian correlation time 

Let us present general analysis of Lagrangian statis- 
tics assuming that the distances between particles are 
smaller than r v . Consider two Lagrangian trajectories 
q(t, ri) and q(t, r 2 ), satisfying the equations d t q(t, j-j) — 
v(t, q(t, j-j)) and the initial conditions g(0, r^) = n. The 
distance between the two trajectories R = q(t,rx) — 
q(t, r 2 ) satisfies the equation d t R — v(t,q(t,ri)) — 
v(t,q(t,r2)) ~ crR, when R is much smaller than the 
viscous length and the velocity difference can be approx- 
imated by the first term in its Taylor expansion. The 
rate-of-strain matrix, 



cr a p(t) = 



dv a (t,r) 



dr 6 



(2.2) 



r=q(t,r 2 ) 



is what determines the deformation of a fluid blob of a 
small size. If the velocity statistics is spatially homoge- 
neous, then the uniform sweeping is statistically irrele- 
vant and only the deformation part of the Lagrangian 
transformation is significant. It is given by R — WRo, 
where the evolution matrix W satisfies 

d t W = aW, W\ t=0 = l. 

The Lagrangian evolution may be thus considered as a 
linear mapping given by the affine transformation W, 
whose statistics is determined by the statistics of a. 
Such Lagrangian mapping can be described universally 
at times much larger than the correlation time, r, of a. In 
the one-dimensional case one finds In W = p = J Q dt'a(t), 
so that at t > t one deals with a sum of large number 
of independent random variables. The probability distri- 
bution function (PDF) of p is completely determined by 
the entropy function S (see e.g. Ref. |5|) 



V(t,p) 



Z(t) 



exp 



-tS 



p-Xt 
t 



(2.3) 



where A is the average value (a) and l/Z is the nor- 
malization factor. At p = Xt the PDF has a sharp 
maximum, described by the central limit theorem, where 
S(x) « x 2 /(2C). Here C is the dispersion of a: 



C 



dt'((a(t)a(t'))) 



This analysis can be generalized to higher dimensions. 
The idea is briefly presented here following the recent 
exposition in ||. At large times, the Lagrangian trans- 
formation can be represented as a stretching or contrac- 
tion along fixed orthogonal directions followed by a rota- 
tion. Indeed, one can represent W as the product MAN, 
where M and N are orthogonal matrices, and A is a 
diagonal matrix |?]J|. At large times the matrix N be- 
comes asymptotically constant, as follows from the Os- 
eledec theorem for W t W 0. Excluding the constant 



matrix N by the proper choice of initial basis one finds 
W = MA. This representation shows that in the frame 
rotating with the fluid blob the Lagrangian mapping is 
just a stretching along fixed directions. Stretching ac- 
cumulates so that the characteristic time of change of 
the matrix A is t while the matrix M changes on the 
much shorter time-scale which is of order of the inverse 
elements of a. Since the time scales of M and A are 
widely different, the matrices are statistically indepen- 
dent. Indeed, changing a at the last stage of evolution 
with duration r one changes M completely, whereas A is 
changed by the amount of the order r/t < 1. It means 
that fixing the value of A(t) does not change the distri- 
bution of M (t) . The matrix M is uniformly distributed 
over the rotation group. The PDF of the eigen values of 
the matrix A = diag [exp(pi), . . . , exp(p c ;)], is given by 



P(t,pO = ie*p [-^(^- Al< 
x0(px - p 2 )--0{pd-i - Pd) ■ 



Pd - x d t 

t 



(2.4) 



Here 9 is the step function that orders the eigen values, 
so that p\ > p2 > •• > pd- The constants A^ (ordered in 
the same way) are the Lyapunov exponents of the flow. 
In principle, they can be expressed via the statistics of a 
(see e.g. 1^,0). We will assume that the Lyapunov spec- 
trum is non-degenerate, i.e. A i > X2 > . . . > X d . The 
normalization factor Z in Eq. ( ^.4[ ) is a function of time. 
Equation (|]4|) shows that at times much larger than the 
correlation time of <r, independently of the statistics of a, 
the statistics of W is characterized by a single function 
S of d variables. This entropy function is convex and 
positive, and it has the expansion 



00 f A A OC q 



(2.5) 



Here Cy is 



near its minimum at X\ = . . . = Xd = 
the covariance matrix of a. Note that we assume the en- 
tropy function to be nonzero at least in some interval of 
Pi which physically means that the flow is random. We 
also assume that the flow is compressible so that S is a 
regular function of all variables (the distribution in the 
incompressible case is obtained by a limiting operation). 
The homogeneous dependence of the PDF on pi and t is 
often sufficient to establish universal statistical laws in- 
dependently of the details of velocity statistics ||[l0 11 1 . 

At large times, V(t,pi) has a sharp maximum at 
Pi = Xit. The existence of generally non-zero Lyapunov 
exponents manifests itself in various growth rates. For 
example, the average logarithmic rate of separation of 
two particles located within the viscous scale of the flow 
is given by Ai, the corresponding growth rate of the fluid 
volume is J2t=i e * c - 

For an incompressible random flow Ai > p2 ], while 
the average value of a is zero (a) — 0. The appearance 
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of non-zero Ai is related to the interplay between rota- 
tional and stretching degrees of freedom [^3|J^Jl^] (for a 
scalar equation one would have Ai = (a) = 0). The 
simplest way to appreciate the existence of positive Lya- 
punov exponent is to consider (following |l^]) an example 
of a saddle-point 2d flow v x — Xx,v y = —Ay where the 
stretching directions satisfy cos</> > (1 + A 2 ) -1 / 2 that is 
their measure is larger than 1/2. 

Compressibility introduces another mechanism of cor- 
relations affecting the stretching. There are more La- 
grangian particles in the contracting regions with V • v < 
0, leading to to the appearance of negative average gra- 
dients in the Lagrangian frame. By isotropy one has 
d{<r a p) = 8 a p(\7 ■ v(t, r)\ r=q{t ro) ). Averaging over the 
volume, one obtains 

(V-«(t,r)| r=g(t)ro) > = J ^y-V-v(t,r)\ r=q{tjro) 
= JyV-v(t,r)(det\ dq{t > ro) 



V • v(t,r) det 



dr 
dq(t,r ) 



q(t,r Q )=r 



dr 



(2.6) 



We observe that this Lagrangian average generally coin- 
cides with the Eulerian average / drV ■ v(t,r)/V only 
in the incompressible case, where it is zero. In the com- 
pressible flow, the integral (|]^) is zero at zero time (when 
we set the initial conditions for the Lagrangian trajecto- 
ries so the measure is uniform back then). The average 
is getting time-independent and negative at times larger 
than the velocity correlation time when the compressed 
regions with negative V • v acquire higher weight than 
the expanded ones. Let us illustrate this conclusion by 
considering the physically interesting case of V • v short- 
correlated in time (to be discussed in much detail below). 
Taking t in (^1]) larger than the correlation time of V • v 
yet small enough to allow for the expansion 



one finds 



1 


dq(t,r) 






dr 





y ■ v{^r)\ r=q{t ro 

1 



2 Jo 



dt'(V -v(t,r)W -v(t',r)) 



Negative (tra) = V • v leads to a suppression of the 
stretching by the velocity field. If one decomposes a into 
"incompressible and compressible parts" 

a aj3 = (^a aj3 - i<5 Q/3 tra^j + i<5 Q/3 tra , 

then from the explicit expressions for Lyapunov expo- 
nents m?]] it is easy to find that the Lyapunov expo- 
nents of the incompressible process a a p — 8 a ptra / d (with 



Ai > 0) get uniformly shifted down by (tra) /d. At a suf- 
ficient degree of compressibility, all the exponents may 
become negative (in the one-dimensional case where com- 
pressibility is maximal one can prove that A < 0, see 
below) . 

Let us stress the difference in Eulerian and Lagrangian 
averages appearing in the compressible case. An Eulerian 
average is uniform over the space while in a Lagrangian 
average every trajectory comes with its own weight deter- 
mined by the local rate of volume change. This difference 
is of an utmost importance in the discussion of backward 
in time Lagrangian statistics to which we pass now. 

We have seen that to find the concentration we must 
consider the evolution of a fluid blob backwards rather 
than forward in time. This is a general situation: to find 
the value of an advected field at the given space-time 
point, r, T, one should consider the Lagrangian trajec- 
tory q(t\T, r) fixed by its final (rather than initial) posi- 
tion: 

d t q(t\T, r) = v(t, q(t\T, r)), q(T\T, r) = r. 

The initial point q(0\T,r) depends on the velocity real- 
ization, so that it is random and not fixed as in the above 
analysis. We denote the Lagrangian quantities related to 
the trajectories fixed by their destination by the tilde 
sign. The strain matrix is defined as 



dv a (t,r) 



drp 



(2.7) 



r=g(t|7> ) 



This must be compared with the definition of a, where 
the initial condition fixes the Lagrangian trajectory. The 
matrices a and a generally have different statistical prop- 
erties. The evolution matrix W is now defined by 

d t W(t\T, r) = aW(t\T, r), W(0\T, r) = 1 

Its value at t — T determines the deformation of fluid 
blobs coming to a fixed final point. Let us relate W to 
W which we now write with a spatial argument 



Wij{t\t',r) = 



dgi(t\t',r) 
drj 



The expression for W(t\T,r) in terms of W(t\t',r) 
is given by W(t\T,r) = W(t\T, r)!^" 1 (0\T, r), so that 
W(T) = W-^T,^). Differentiating the identity 
q(T\0,q(0\T,r)) = r one finds W(T\0, q(0\T, r)) = 
W~ l {0\T,r), which relates W and W: 

W(T\T,r) = W(T\0,q(0\T,r)). 

To express the statistics of W in terms of the La- 
grangian characteristics we use the same transformation 
that we used for transforming the average of V • v in the 
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Lagrangian frame to the usual Eulerian average. Namely, 
for the average over the volume of the flow of any function 
/ one has 

(f{W(T)})= j ±lf {W (T\0, q (0\T,ro))} 



dx 



f{W(T\0,x)}det 



dq(T\0,x) 



dx 



= (f{W(T\0,x)}detW(T\0,x)) . 



(2.8) 



Again, Lagrangian and Eulerian averages coincide for in- 
compressible flow, when detW = 1. In general, it is 
necessary to account for the local volume change when 
passing from one average to another. Since detiy = 
exp^/?i, then considering passive fields one has to 
take the following probability distribution of stretch- 
ing/contraction eigen values 



Pi - Aii pd - X d t 



f 



t 



x6»(pi - p 2 ) ■ ..0{p d -i 







(2.9) 



Here S and Z are the same as in J2.4|) that is what one 
can measure in studying (forward-in-time) particle dis- 
persion. Note that the correct normalization of V is guar- 
antee d by the volume conservation (det W) = 1 following 
from (U) with / = 1. 

The Lyapunov exponents Aj of W are determined by 
the extremum of the exponent: 



A; 



Aj + Vi 



(2.10) 



where are determined from dS{y\, ..,y d ) / dyi = 1. 
An important remark is that yi cannot generally be ex- 
pressed via the Lyapunov exponents A; only. They de- 
pend on the form of the entropy function S and hence 
on the details of velocity statistics. Indeed, every trajec- 
tory comes with its own weight determined by the local 
rate of volume change. The consequence is that pas- 
sive fields behavior in a compressible flow does not enjoy 
the same degree of universality as in the incompressible 
case (when, for instance, the growth rate of the magnetic 
fluctuations is determined solely by the spectrum of the 
Lyapunov exponents Xi and is independent of the form 
of the entropy function [fl0| ). 

Even though we will use only the properties of the 
matrix W, we also mention the probability distribution 
function of the eigen values of matrix VK(0|T, r) which 
directly determines the evolution backwards in time 



i=i \ 



Pi - A d f 
t 



-Pd - Mt 
t 



0{pi - P2) ■■■ 0{p d -i - Pd) 



which follows from the above results using W(0|T, r) = 
W (T\T, r). It follows that the backward-in-time Lya- 
punov exponents are given by —Xi and not by the naive 
guess — Aj, which holds only in the incompressible case. 
In particular, particles diverge backwards in time with 
exponent — A^. 

The difference between Lyapunov exponents Xi and Xi 
can be illustrated in the one-dimensional case where their 
signs are definite and opposite. Indeed, the conservation 
of the total fluid volume together with spatial homogene- 
ity imply 



p-tS 



p-Xt 
t 



= 1. 



(2.11) 



Let us consider this identity at large times when the 
integral can be calculated by the saddle-point method. 
First, we note that at large t the normalization factor 
Z = J dpenp(-tS) is determined by the region, where 
S(x) cx x 2 so that Z cx i 1 / 2 . Next, one can rewrite the 
integral ( |2.1l| ) as tZ^ 1 J dx exp[i(A + x - S(x))]. It is 
determined by the point where x — S{x) is maximal. 
Since x — S{x) takes positive values near x — we con- 
clude that x*— 5(x») > 0. Therefore A = S(x*) — x* < 0, 



for the integral ( p. 11 ) to be time-independent. 

On the contrary, A is positive. Defining p + 
S((p — Xt)/t) = S((p — Ai)/t) one has the condition 
J dpexp[—p — tS(p/t — A)] = 1 which gives A > 0. 

The above results are generalized to higher dimensions 
as follows. Since we consider the flow to be contained in a 
fixed volume then (det = 1, and in the same manner 
one finds that the mean logarithmic rate-of-change of the 
volume elements Aj < (which, in particular, implies 
A3 < 0). From the explicit expressions for Xi one 
finds — ( V • v) , so that we have proved that the La- 
grangian average (V • v) is nonpositive, the result stated 
above on physical grounds. The corresponding inequality 
on Ai is ^ Aj > (implying Ai > 0). We arrive at a some- 
what surprising conclusion that for the Lagrangian dy- 
namics one has average compression of volumes, whereas 
passive fields rather feel average expansion. The physi- 
cal meaning of this effect is transparent: as we go away 
(either forward in calculating Xi or backwards in calcu- 
lating Aj) from the moment where we imposed a uniform 
Lagrangian measure, the volume rate-of change is getting 
negative in a fluctuating compressible flow. To avoid mis- 
understanding, let us stress that for a physical quantity 
x(t) (volume of a fluid element in this case) the conser- 
vation of the mean value (x(t)) does not contradict to a 
nonzero rate of change t~ l (lu.x(t)}. 

The general considerations can be illustrated using a 
particular case of the velocity statistics, the Kraichnan 
model of a short-correlated Gaussian velocity with the 
variance 

(v a (t, r)v (0, 0)) = [V 5 af3 - JC af3 {r)] 5{t) , (2.12) 
Zap = D [(d + 1 - 2T)S a pr 2 + 2{dT ~ l)r a r p ] . (2.13) 
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respectively, it is thus the degree of compressibility that 
may vary between and 1. The quadratic dependence of 
the correlation function on the coordinate corresponds to 
the expansion of the velocity difference we made above. 
For a velocity defined by Eqs. ( |2.12 ) and ( [2.13] ), a 
straightforward calculation gives 



Xi/D = d(d + l- 2i) - 2T[d + (d- 2)i\ 



(2-14) 



In the incompressible case, T = 0, this formula has been 
derived in fl4|| . For a general compressible case, Ai has 
been derived in (i) , where it has been also observed that 
Ai chang es si gn at F = d/A. The entropy function has 
the form (2J3) for arbitrary values of x. One can also find 



C i3 = AD{[d + T{d - 2)]5ij - 1 + Td} 



We see from ( |2.14| ) that compressibility indeed diminishes 
the Lyapunov exponents. It is interesting to compare 
this with the influence of compressibility on Lagrangian 
dynamics in a multiscale velocity: there, the Lagrangian 
trajectories either explosively separate or implosively col- 
lapse depending on whether the degree of compressibility 
is small or large respectively [ [l5f . 

The Lyapunov exponents Ai that govern the behavior 
of the passive fields are enhan ced by compressibility since 
Di are positive (see Eq. ( 2.1 0|) ) . For the Kraichnan model 
one has yi — J^j Cij/2, so that 

Xi/D = d(d + l- 2i) + 2T[d 2 -(d- 2)i - 2} . (2.15) 



For d = 2, 3 one has 

Ai = 2L>(1 + 2r) , A 2 = -2D{1 - 2T) , 
Ai = 6L»(1 + 2r) , A 2 = 10LT , A 3 / = 



-2L>(3 - 4T) . 



The compressibility, T, is identically equal to unity in 
d = 1, where instead of Eq. ( [2.15 ) one should write 
K = Dx 2 . Then, Eqs. @.14| ) and fl2.15|) are replaced by 
A = — D and X — D respectively. 

Comparison of (2.14) and ( |2.15 ) shows that Ai = 
— Xd+i-i- This relation is due to time-reversibility of the 
short-correlated velocity. In particular, Ai and A^ change 
sign at the same degree of compressibility T — d/4. This 
is peculiar for a short-correlated case and does not hold 
for an arbitrary velocity statistics. In other words, the 
change of the regime from stretching to contraction in the 
forward Lagrangian dynamics does not generally corre- 
spond to the change of the regime in the passive fields, 
which are related to the backwards in time Lagrangian 
dynamics. 



B. Applicability of the ideal case approximation 

Let us now determine the domain of validity of the 
ideal case approximation. Our starting point will be the 



dynamic expression for the Green's function that can be 
derived explicitly in this case. To find when the concen- 
tration starts to be determined by the particles that were 
initially separated by a distance larger than the viscous 
scale, we must analyze the support of the Green's func- 
tion, G(t,r\t' = 0,r'), as a function of the initial coordi- 
nate, r' . At k — the dynamics is purely Lagrangian, so 
that 



G(t,r\0,r') = 



1 



det W(t, r) 



S(r'-q(0\t,r)). 



This formula has a clear meaning: the change of the con- 
centration at a point is completely determined by the 
volume compression factor along the Lagrangian trajec- 
tory. Small diffusion is equivalent to adding a Brownian 
motion to the velocity. It leads to a smearing of the re- 
gion around the Lagrangian trajectory from which par- 
ticles come to the observation point. As a function of 
the initial arguments, the Green's function satisfies the 
Hermite-conjugate evolution equation 



d t ,G +(v(t', r 1 ), Vr/)G = -kV^G. 



(2.16) 



As long as the support of G is much smaller than the 
viscous scale, one can expand the velocity in the vicinity 
of the Lagrangian trajectory in the Taylor series. The 
homogeneous component is excluded by passing to the 
moving frame and the first non-trivial term contains a: 



d t >G + & afJ {t'y fJ V a G = - K V 2 r ,G . 

This equation can be solved in the Fourier space 
1 



(2.17) 



G(t,r|0,r') = 



detW{t,r) 



dk 



■ exp 



ik-[r' -q(0\t,r)]- 



k l Ik 



(2.18) 



I = 2k dt' W- l {lI\t,r)W- ht {t'\t,r) . 
Jo 

The matrix / is the inertia tensor of a patch of particles, 
evaluated at t = 0, provided the patch is a sphere with 
the center at the point r at time t. The particles perform 
independent Brownian motions together with the La- 
grangian motion in the same velocity field (cf. 0). Let us 



stress that no averaging has been performed in Eq. (2.1J ) 
and therefore the expression for G is purely dynamical. 
We observe that the size of the region which makes the 
main contribution to the concentration at a point grows 
as the largest eigen value of the matrix I, i.e. the square 
of the linear size grows as k L dt' exp[— 2pd{t')]. Since the 
diffusionless consideration is valid as long as the largest 
size of the rd-volume is smaller than the viscous scale, 



the applicability condition of Eq. (2. IS) is 



G 



k I dt' eiq>[-2p d (t')] « 
Jo 



(2.19) 



Below, we will refer to the configurations of velocity 
with decreasing pd(t) as the contracting configurations. 
Indeed, for such configurations, particles at the observa- 
tion point are brought together from larger regions. On 
the other hand, for the configurations with increasing 
p~d(t), the concentration is determined by the particles 
initially belonging to the region of the size of the or- 
der yj k/X. Such configurations can be called diverging 
because the particles in the vicinity of q(0\t, r) diverge 
exponentially. 

Formula ( ]2.18| ) gives n(t,r) = J G(t,r\0,r')dr' = 
1/ det W(t, r). This is exactly the same expression as 
for non-diffusing particles. For the moments of the con- 
centration one finds 



J d Pi exp[-(a- 1) 



Pi]P(t, Pi ) 



(2.20) 



with P(t,pi) given by Eq. ( |2.4| ). Note that the growth 
function 7l(«) in the Lagrangian frame is obtained by 
a mere shifting of the argument of the Eulerian growth 
function: 



(n a (q(t,r),t))=(n a+1 (r,t)) 



(2.21) 



Integral Q2.20| ) can be calculated using the saddle-point 
approximation. The saddle-point value of pd given by 
—Cat, where c a is an a— dependent constant, i.e. the 
condition of applicability is k J q dt 1 exp[2c a t'] <C r\. Us- 
ing convexity of the entropy one can show that large neg- 
ative moments have negative c a and therefore the above 
condition becomes time- independent. This can be sim- 
ply seen noting that the averaged quantity exp[— a^2 pi] 
favors positive pi at negative a. Therefore, the diffusion- 
less result is always correct for large negative moments. 
Diffusion cannot stop the formation of void regions with 
few particles inside. 

On the other hand, from expression J2.2C| ) one can see 
that for a > any growing moment of n must be deter- 
mined by a positive c a , otherwise Pi > Pd > 0- Gen- 
erally, one can assert the existence of the boundary ah, 
such that — oo < ccb < 1- For a < ah, the saddle-point c a 
is negative and the corresponding moment behaves as in 
the diffusionless case for all times, whereas for a > ah the 
diffusionless approximation breaks down at large times. 
Since the moments with a < ah arc determined by 
the configurations on which the rd-volume is compressed 
backwards in time, one expects that ah is a monoton- 
ically increasing function of the velocity compressibility 
(as measured by ^ Aj). For example, in the framework of 
the Kraichnan model one has «b = (r — 4,d)/(2T(d + 2)). 
We will refer to the ah < case as the weakly compress- 
ible case, and < ah < 1 as the strongly compressible 
one. It can be verified that this corresponds to the cases 



of Xd < and > respectively. The same is valid for 
an arbitrary velocity statistics. 

In fact at any time t one can consider the contribution 
of diverging configurations, so to say, the "ideal fluid con- 
tribution" : 



/id = 



Pd>0 



dpi exp —(a -I) 22 Pi 



ts 



(2.22) 



Since the smallest size cannot be smaller than rj, it is 
necessary to introduce here the c utoff at pd = 0. It is 
clear that (n a ) > (n a )id- Integral ( 2. 22] ) has exponential 
time-dependence, (n a )id oc exp[7;d(o:)t]. Note that due 
to the constraint pd > the growth function 7;d is differ- 
ent from 7(a). For a < ah the saddle-point is inside the 
domain of integration at all times. On the contrary, for 
a > ah, integral (2.22) is determined by the boundary, 
Pd = 0. 

In the weakly compressible situation, ah < 0, one has 
7id(c>!b) = 7(c*t>) > 0. From the continuity of 7i d (a) we 
conclude that 7i d (a) > for a < a' h < 0, where a' h is de- 
fined as 7id(a4) = 0- The inequality a' h < follows from 
7id(a) < 7(a) and 7(0) = 0. Therefore, the moments of 
the order a < a' h < satisfy (n a ) > exp(j' a t) with posi- 
tive "f' a at all times (in fact, asymptotically the equality 
holds as mixing configurations can only lead to a growth 
slower than exponential, see below). It means that these 
moments become infinite in the steady state, which cor- 
responds to the formation of the power-law asymptotic 
behavior for the PDF of the concentration near n = 0: 
7 5 (n) cx n _ "t> _1 . Diffusion does modify the growth of the 
moments with ah < a < a' h , but the time dependence 
remains exponential. 

In the strongly compressible case, ah > 0, one can im- 
mediately conclude that the moments with a < grow 
exponentially with the ideal fluid exponents. For a > 
the inequality (n a ) > (n Q )id leads to no interesting con- 
clusions at large times. The asymptotic behavior of the 
concentration PDF at n — > is V(n) oc nT 1 . 

Let us now consider the moments to which the main 
contribution is made by the contracting configurations 
(i.e. the moments with a > 1). The cutoff time t* 
for the ideal growth is determined from the condition 
exp(2c Q i*) ~ (r v /r d ) 2 , which gives t* a = c^ 1 ln(r v /V d ). 
This expression is exact in the limit of large Schmidt 
numbers. Note that the cutoff time depends on the or- 
der of the moment. For a quadratic in a entropy (e.g. 
for the Kraichnan model), c a is a linear function of a. 
The steady-state dependence of the moments on Schmidt 
number can be estimated from below by (r v /rd) _7 ^ Q ^ c ° . 
One can expect that the dependence of the steady-state 
moments on the Schmidt number is linear for large a, 
since 7 Q ~ ac a . This can be seen from the saddle-point 
expression for the moment. The linear dependence signi- 
fies less intermittent tail of the PDF as compared to the 
evolution problem. 
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III. SATURATION OF GROWTH DUE TO 
DIFFUSION 



To analyze the behavior of the moments at larger times 
one should distinguish two cases: the Re ~ 1 case when 
the velocity correlation length L is of the order of r v , 
and the case of Re 3> 1, when L ^> r v . In the first case, 
advection becomes equivalent to the usual diffusion at 
the scales larger than r v . The moments get saturated 
and are given by the corresponding power of r v /rd. In 
the large Re case the velocity divergence is correlated 
at scales much larger than r v . The correlation between 
different viscous domains decays as a power of the dis- 
tance between them. The configurations that coherently 
bring together different viscous domains determine the 
moments at this stage of evolution. 

The moments continue to grow (in a power-law fash- 
ion) only for a particular case when the compressible cor- 
rection to the velocity has the same scaling as the incom- 
pressible velocity. Since the compressible part is propor- 
tional to (uV)u, then u and v have the same scaling only 
for a smooth velocity, Su oc r, which has been studied in 
Sec. [n]. Note that up to logarithmic corrections this is 
true for a vorticity 2d cascade as well. However, for the 
turbulent velocity in the energy cascade, the velocity u is 
non-smooth, hence the compressible part has a different 
scaling. For example, the Kolmogorov phenomenology 
gives Su cc r 1 / 3 ), so that Sv oc r -1 / 3 . It means that 
the compressibility is most important at small (viscous) 
scale so that the growth has to saturate and the level of 
fluctuations should not depend on the Reynolds number. 

Unfortunately we still lack the formalism to describe 
Lagrangian statistics in the inertial interval with the 
same degree of universality as in the viscous interval. 
Nevertheless, to understand the most essential proper- 
ties of the concentration fluctuations, one can use the 
simplest velocity statistics. We assume that the velocity 
is statistically isotropic, Gaussian, and has zero corre- 
lation time p6| . The pair-correlation function is given 
by 



(v a (t, r, )^(0, 0)) = 8(t)[V 5 a p - K a p{r)\ , 



(d - I)/C a/ 3 



( r d+i u y 



r 2 S a f 



(r 2 u)' 



(3.1) 



Note that c = for incompressible flows. 

We will assume that u and c have a regular expansion 
at r « r v , that is u(r) w it(0) + u"(0)r 2 /2 + . . ., and 
c ps c(0) + c"(0)r 2 /2 + . . .. In the intermediate region 
r v <^ r <^ L, the functions u and c behave in power- 



law manners. However, due to relation (1.1), the scaling 
exponents of u and c are generally different. At large 
scales, r ^> L, the function u scales as r~ 2 . To guarantee 
that the variance of the velocity is positive, c must vanish 
faster than r~ 2 at r ^> L. 



A convenient measure of compressibility is given by the 
ratio e = c/u. We will denote by eo the value of e for the 
smooth velocity (u,c constants). Despite all these crude 
simplifications the model enables to see the most inter- 
esting features of the growth and is used to illustrate the 
conclusions we believe to be model- independent. 



A. The two-point correlation function 

In this subsection we study the two-point correlation 
function of the concentr ation , /(r) = (n(O)n(r)). In the 
framework of the model (3T) it satisfies the closed equa- 
tion 



d t f = V a V /3 (Kafif) + 2kV 2 /. 



(3.2) 



Let us first show that the dynamics of / is relaxational 
and then find the stationary solution to which it con- 
verges. For this purpose we must analyze the spectrum 
of the differential operator on the right-hand side of Eq. 
( |3.2[ ). The eigen functions must be regular at small dis- 
tances. The boundary condition at large distances fol- 
lows from the fact that the correlation function tends to a 
constant, equal to 1 (?1q in the dimensional units). There- 
fore one must require that the eigen functions do not grow 
at infinity. The operator has the form of a d— dimensional 
Fokker-Planck operator, so that one can expect that it 
has no positive eigen values. Due to spherical symmetry 
of / on can rewrite Eq. (13.21) in the spherical coordinates 



dtf = r 1 - d £ FP 



„d-l 



/) 



Cfp = d r [r 2 u + 2k] e *<9, 



— rexp 



rc(r) dr 
r 2 u(r) + 2k 



(3.3) 
(3.4) 

(3.5) 



Note that C fp has the form of a one-dimensional Fokker- 
Planck operator. Now we can show that the operator C 
on the right-hand side of Eq. (3.2) has no positive eigen 
values, i.e. all the eigen functions of the operator sat- 
isfy £Je = —EJe with E > 0. Indeed, the evolution 
operator becomes proportional to a Laplacian at r ^> L. 
Therefore the negative energy eigen functions have ex- 
ponential behavior at infinity. The boundary condition 
at infinity ensures that only exponentially decaying so- 
lutions are allowed. To show that for such solutions the 
boundary condition at r — cannot be satisfied, we write 
the identity 



dr f E e*r^£f E = -E dr f E e*r^ 
i« Jo 

'dr (r 2 u + 2k) e~* (d r e^r d ' 1 f E f < 



which proves that E > 0. We used integration by parts, 
which is possible only for functions decaying at infinity. 
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Let us now show that in fact the spectrum is contin- 
uous, covering the interval [0, oo). For definiteness we 
assume that < r v < L. Let us consider the equation 
at r < r v where u = u(0) and c = eit(0): 



(l + x 2 )/" + 



(d+l + e)x + 



1 



/' + (de + A)/ = 0. 



Above we have assumed that the distance is measured in 
units rj, where r 2 = 2k/u(0). We have also introduced 
A = E/u(0). The solution of this equation satisfying the 
correct boundary conditions is 




where /i = y/(d — e) 2 — 4A. Next, let us consider the in- 
ertial interval. Considering for simplicity u = Dr^ 1 and 



c = e y Dr 7 we find 



r 2 f" + (d+l- 7 + e)r/' + 
The solution of this equation is 

/ = C 1 {E)r-^ d+e -^l 2 J u 
+C 2 (E)r- [ - d +^V 2 N iy 



(d-7)e + 



Sr 7 



/ = 0, 



'2\/i? 



.7/2 



Wd 1 



.7/2 



(3.6) 



where v = (d — 7 — e)/^. Finally, at r ^> L one can 
use the function ( p.q ) substituting 7 = 2, e = and Vq 
instead of D 

f = C 3 (E)r 1 - d / 2 J d/2 _ 1 { x l§r 




+C i {E)r 1 - d / 2 N d/2 _ l 1 r 

We observe that unlike the case of E < 0, for E > 
there is no additional restriction on the coefficients of 
the eigen function. Thus the matching problem can al- 
ways be solved. The matching at r ~ r v fixes the ratio of 
constants C\/C 2 , and the matching at r ~ L v fixes the 
ratio C3/C4. We conclude that the spectrum is positive 
continuous and non-degenerate. In fact, this property 
holds for any relation between and r v . 

Thus at large times / must converge to /o = / s t- Note 
that decay at large times is prohibited by the inequality 
f(t, 0) = (n 2 ) > (n) 2 = 1. Let us now find the stationary 
solution which satisfies 



d r ([r 2 u + 2k] e~*d r (eV-7 st )) = 0. 



(3.7) 



The function must approach the square of the average 
concentration at infinity and be regular at zero. It is 



easy to see that the solution satisfying these conditions 
is the "zero flux" solution, proportional to exp[— $]r 1_rf . 
It is given by 



/si 



exp 



r'c{r')dr' 



r' 2 u{r') + 2k 



(3.8) 



The integral in Eq. (3.8) converges, because the function 
c decays faster than r~ 2 at r ^> L. Since r 2 u ~ Vq at 
r>L, we obtain the following asymptotic expression at 
these scales 



/si 



I \ [ - 

r 



where we assumed that c decays as r~ 2 ~ a , a > 0. 

The behavior of / s t in the inertial interval, r v <C r <C 
i, crucially depends on whether u and c have the same 
scaling exponents. If they do, f s t(r) behaves as a neg- 
ative power of r. The single point correlation function, 
(n 2 ), which can be estimated as / s t(fd), is then propor- 
tional to a positive power of the Reynolds number. If, 
however, the velocity is not smooth in the inertial inter- 
val, u and c have different scaling exponents, and Eq. 
( |3.8[ ) shows that the main growth of / st occurs below the 
viscous scale. Hence, (n 2 ) is independent of the Reynolds 
number. For example, for the Kolmogorov scaling, the 

8/ 3/3 - 



solution (3 
r <C r v . 



has the form ln/ st oc a 



-4/3, 



at 



Therefore, fluctuations of the concentration are mainly 



produced in the interval of scales r 



< 



where the fluid 



velocity is smooth (i.e. in the viscous interval or in 2d 
vorticity cascade). One can write estimates 



(3.9) 



Since is by definition larger than a, significant fluc- 
tuations are possible only for very heavy particles with 
13 m 2p/3 Pp < (a/r v ) 2 ln 1/2 (r v /r d ). 

To conclude, the fluctuations of concentration grow 
exponentially in any random compressible flow with a 
nonzero sum of Lyapunov exponents until this growth is 
restricted by finite-size effects (diffusion or discreetness 
of the particles). It is interesting to note that n can 
be considered as the density of the fluid itself, so that 
the finite-size effects are absent. Since density pertur- 
bations does not grow unlimited, one concludes that the 
phenomenon described here can take place only as a tran- 
sient process when, for instance, large-Mach random flow 
with almost homogeneous density was initially created. 
In a stationary turbulence of the compressible fluid, the 
back reaction of density fluctuations on the flow stop the 
growth of the density fluctuations (we are indebted to V. 
Lebedev for this remark). 
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